#include "cells.hpp"
#include "mesh.hpp"
#include "geometry.hpp"
#include "topology.hpp"
#include "signalCodes.hpp"
#include "exceptions.hpp"
#include "util.hpp"
#include <utility>
#include <cmath>
#include <cassert>

namespace hms
{

Cells::Cells( const Mesh& mesh ) :
	CellsBase<Cells>(mesh)
{}


auto Cells::centroidsImp() const -> CentroidsReturnType {
	return m_centroids;
}
auto Cells::areasImp() const -> AreasReturnType {
	return m_areas;
}

auto Cells::circumfsImp() const -> CircumfsReturnType {
	return m_circumfs;
}

auto Cells::nDomainCellEdges() const -> const ArrayX<char>& {
	return m_nDomainCellEdges;
}
auto Cells::domainCellEdges() const -> const Array4Xi& {
	return m_domainCellEdges;
}
void Cells::nDomainCellEdges( ArrayX<char>&& nEdges ){
	assert( m_nDomainCellEdges.size() == 0 );
	m_nDomainCellEdges = std::move(nEdges);
}
void Cells::domainCellEdges( Array4Xi&& edgeIndices ){
	assert( m_domainCellEdges.size() == 0 );
	m_domainCellEdges = std::move(edgeIndices);
}
void Cells::centroids( Matrix2Xs&& c ){
	assert( m_centroids.size() == 0 );
	m_centroids = std::move(c);
}
void Cells::areas( ArrayXs&& a ){
	assert( m_areas.size() == 0 );
	m_areas = std::move(a);
}

void Cells::storeData(
	ArrayX<char>&& nEdges,
	Array4Xi&& edges,
	const Ref<const Array4Xi>& domainCellVertices
){
	nDomainCellEdges( std::move(nEdges) );
	domainCellEdges( std::move( edges) );
	storeAreasAndCircumfs(domainCellVertices);
}


auto Cells::nDomainCells() const -> Index {
	return m_domainCellEdges.cols();
}
auto Cells::nGhostCells() const -> Index {
	return mesh().nBoundaryEdges();
}
auto Cells::nCells() const -> Index {
	return nDomainCells() + nGhostCells();
}

#define ASSERT_INDEX \
	assert( cell >= 0 && cell <= nCells() && "Invalid cell index");

auto Cells::area(Index cell) const -> scalar {
	ASSERT_INDEX
	return m_areas[cell];
}

auto Cells::centroid(Index cell) const -> CentroidType {
	ASSERT_INDEX
	return m_centroids.col(cell);
}

auto Cells::nEdges(Index cell) const -> char {
	ASSERT_INDEX
	return isGhost(cell) ? 1 : m_nDomainCellEdges[cell];
}

auto Cells::edge(Index cell) const -> Index {
	ASSERT_INDEX
	assert( isGhost(cell) &&
		"Cells::edge(Index) is for ghost cells only!"
		"For domain cells, use Cells::edges(Index) instead!");
	
	return mesh().nInnerEdges() + cell - nDomainCells();
}

auto Cells::edges(Index cell) const -> EdgeIndices {
	ASSERT_INDEX
	assert( !isGhost(cell) &&
		"Cells::edges(Index) is for domain cells only!"
		"For ghost cells, use Cells::edge(Index) instead!");
		
	return m_domainCellEdges.col(cell).head(m_nDomainCellEdges[cell]);
}

auto Cells::isGhost(Index cell) const -> bool {
	ASSERT_INDEX
	return cell >= nDomainCells();
}

auto Cells::domainCell(Index cell) const -> Cell {
	assert( isGhost(cell) );
	return mesh().cell( mesh().edge( edge(cell) ).domainCell() );
}



auto Cells::isBoundary(Index cell) const -> bool {
	ASSERT_INDEX
	/* return early if it's a ghost cell */
	if (isGhost(cell)) return false;
	/* otherwise look at edges. If one is a boundary edge, return true */
	for ( char i {0}; i<m_nDomainCellEdges[cell]; ++i )
		if ( mesh().edge( m_domainCellEdges(i,cell) ).isBoundary() )
			return true;
	return false;
}

auto Cells::circumf(Index cell) const -> scalar {
	ASSERT_INDEX
	return m_circumfs[cell];
	// Index domainCell;
	// if ( !isGhost(cell) )
	// 	domainCell = cell;
	// else
	// 	domainCell = mesh().edge( edge(cell) ).domainCell();
	
	// return mesh().edges().lengths()( edges(domainCell) ).sum();
}

/**
 * @brief Collects neighbouring cell indices by looking up parent cells of
 * constituting edges. Collection is performed when function is called!
 * @param cell The cell index.
 * @return an ArrayXi whose length is equal to the number of neighbours.
 * Ghost cells have only one neighbour, in which case the return array
 * has only one entry.
 */
auto Cells::neighbours(Index cell) const -> ArrayXi {
	ASSERT_INDEX
	ArrayXi neighbours ( this->nEdges(cell) );
	if ( cell < nDomainCells() ){
		for ( char edge {0}; edge < m_nDomainCellEdges[cell]; ++edge ){
			Array2Xi::ConstColXpr parents {
				mesh().edge( // edge index
					m_domainCellEdges.col(cell)[edge]
				).cells() // parent cells
			};
			neighbours( edge ) = parents[0] != cell ? parents[0] : parents[1];
		}
	} else
		neighbours(0) = mesh().edge( this->edge(cell) ).cells()[0];
	return neighbours;
}

/**
 * @brief The Cells object does not store its vertices. This routine
 * uses the stored edge indices instead, to fetch the vertex indices through
 * the parent mesh. Therefore, this routine has a certain cost.
 * @param cell Index of the cell
 * @param ordered flag for whether or not the output shall be ordered counter-
 * clockwise
 * @return a std::vector<Index> containing the vertex indices of the cell.
 * Can be ordered counter-clockwise, if the respective parameter is set to true.
 */
auto Cells::vertices(Index cell, bool ordered) const -> ArrayXi {
	ASSERT_INDEX
	ArrayXi verts;
	if ( isGhost(cell) )
		verts = mesh().edge( this->edge(cell) ).vertices();
	else
		verts = edgeLoopToVertices(
			mesh(),
			cell
		);

	if (ordered)
		orderVertices(
			mesh(),
			verts,
			VertexOrdering::counterclockwise
		);

	return verts;
}



/**
 * @brief Calculates and stores areas and centroids of all cells.
 * 
 * When a Cells object is initialised with edge indices, which 
 * occurs after vertices and edges are constructed as part of the mesh, then
 * the available information can be used to calculate areas and centroids of
 * each cell. This calculation is coordinated by this routine.
 */
void Cells::storeAreasAndCircumfs(
	const Ref<const Array4Xi>& domainCellVertices
){
	/* first, resize centroid and area matrices to cell count */
	m_centroids.resize( NoChange, nCells() );
	m_areas.resize( nCells() );
	m_circumfs.resizeLike(m_areas);
	
	/* area and centroid calculation both require vertices without duplicates.
	 * With those, the calculation routines are available from geometry.hpp */
	if ( domainCellVertices.cols() == nDomainCells() ){
		/* if domain cell vertices were passed in, they don't need to be
		 * retrieved and ordered */
		mesh().forEachDomainCellIndex( [&]( Index cell ){
			m_areas[cell] = verticesToArea(
				mesh(), domainCellVertices.col(cell).head( nEdges(cell) )
			);
			m_centroids.col(cell) = verticesToCentroid(
				mesh(), domainCellVertices.col(cell).head( nEdges(cell) ),
				m_areas[cell]
			);
		} );
	} else {
		mesh().forEachDomainCellIndex( [&]( Index cell ){
			ArrayXi verts { this->vertices(cell, true) };
			m_areas[cell] = verticesToArea( mesh(), verts );
			m_centroids.col(cell) = verticesToCentroid(mesh(), verts, m_areas[cell]);
		} );
	}
	mesh().forEachDomainCellIndex( [&]( Index cell ){
		m_circumfs[cell] = mesh().edges().lengths()( edges(cell) ).sum();
	} );

	/* for ghost cells, the area is the same as that of the neighbouring
	 * domain cell, and the centroid is similarly the domain cell's centroid
	 * mirrored along the shared boundary edge. */
	mesh().forEachGhostCellIndex( [&]( Index cell ){

		Edge edge { mesh().edge(
			cell - mesh().nDomainCells() + mesh().nInnerEdges()
		) };

		#ifdef MESH_SQUARE_GHOSTS
		/**********************************************************************/
		/* option square ghost cells */
		/**********************************************************************/
		Matrix2s centroids { edgeToSquareCenters( mesh(), edge.index() ) };
		/* For creating a square from one edge, there are two options in.
		 * Let the vertices be A and B and the edge be AB, then the two
		 * remaining corners to form a square can either lie both on the left 
		 * or both on the right of AB. Thus there are two possible centroid points
		 * C1 and C2. Let C1 be on the left and C2 be on the right of AB.
		 * We already know the centroid point P of the neighbouring cell to which
		 * the edge originally belonged. P and the new centroid point must be
		 * on opposite sides of AB.
		 * edgeToSquareCenters() returns C1 and C2, with C1 always lying left
		 * and C2 always lying right of AB.
		 */
		m_centroids.col(cell) =
			whichSide(
				mesh(),
				edge.index(),
				m_centroids.col( mesh().edge( edge.index() ).cells()[0] )
			) == Side::right ?
			centroids.col(0) :
			centroids.col(1);
		
		/* the area is simply the square of the edge length */
		m_areas[cell] = pow( edge.length(),2 );
		/* the circumference is then four times that edge length */
		m_circumfs[cell] = 4 * edge.length();

		#else
		/**********************************************************************/
		/* option same area as domain cell, mirrored centroid */
		/**********************************************************************/

		/* simply copy the area and circumference over */
		m_areas[cell] = m_areas[ edge.domainCell() ];
		m_circumfs[cell] = m_circumfs[ edge.domainCell() ];

		/* project the domain cell centroid onto the edge, resulting in the
		 * point on the edge from which a line normal to the edge goes through
		 * the domain cell centroid.
		 * The ghost cell centroid lies on the other side of the edge, on that
		 * same line, in the same distance */
		Matrix2Xs::ColXpr c { m_centroids.col( edge.domainCell() ) };
		m_centroids.col(cell) = mirrorAlongEdge(edge, c);
		
		#endif
	} );
}

} // end namespace hms

#undef ASSERT_INDEX